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ABSTRACT 


An  improved  method  to  control  error  in  the  classical 
Runge-Kutta  numerical  integration  scheme  has  been 
uncovered.  The  method  is  readily  adaptable  to  the 
solutions  of  the  differential  equations  of  notion  of 
an  unguided  rocket.  This  report  presents  the  theo¬ 
retical  aspects  behind  this  method  and  an  evaluation 
of  its  efficiency. 
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INTRODUCTION 


Most  mathematical  models  for  the  simulation  of  the  trajectory  of  a  bal¬ 
listic  rocket  employ  the  Runge-Kutta  numerical  integration  scheme  with 
the  extrapolation-to-zero-grid  error  control  criterion  of  Richardson  [1]  . 
This  criterion  is  based  on  a  comparison  of  the  Results  obtained  by  nu¬ 
merically  integrating  with  step  sizes  of  h  and  .  This  involves  eval¬ 
uating  the  functional  equations  eleven  times  for  a  single  integration 
step,  seven  of  the  evaluations  being  performed  solely  for  the  purpose  of 
error  control . 

An  alternate  technique  for  error  control  was  suggested  by  Earnest  [2]  . 
This  report  investigates  Earnest's  technique  and  compares  the  efficiency 
of  this  technique  to  that  of  Richardson. 


DISCUSSION 


Suppose  one  has  a  system  of  first-order  ordinary  differential  equations 
of  the  form 


dyi  , 
“37  =  yi 


f i  (x,  yx. 


....  yN)  i  =  i, 


,  N  (1) 

The  Runge-Kutta 


with  initial  conditions  y^(Xy)  =  >\(),  i  =  1, 
fourth-order  method  is  an  algorithm  designed  to  approximate  the  Taylor 
series  expansions 


yi (x0  +  h)  =  yi0  *  h  yiO  *  ^ 

4.  . 

of  (I)  and  has  the  general  form 


iO 


y-(xn  ♦  h)  =  y.0  -  aK.n  ♦  bK.,  ♦  c  K.,  ♦  dK. 


iO 


il 


i3 


(2) 


(3) 


where 


Ki0 

Kil 

Ki? 

K.„ 

10 


h f i ( x0 ’  y10' 

•  •  *  1 

=  hf.(x0  ♦  lh, 

yio 

*  "  ho 

=  hf .  (x.,  ♦  nh, 
i  0 

yio 

*  ■>  ho 

•  hf.(x()  ♦  ph, 

y10 

v  ho 

1. 


y*. 


;o 


+  m 


♦  r  K 


IT 


yN0 


♦  q  K, 


NO 


(4) 

♦  r  Si 


♦  s  k}1  ♦  t  K12, 


...) 


For  purposes  of  completeness  a  brief  discussion  of  the  determination  of 
the  above  parameters  is  included  in  Appendix  A,  the  reader  is  referred 
to  [3]  for  a  more  thorougli  discussion. 


THE  ERROR  CONTROL 


1 .  Derivation 

One  assumes  a  truncation  error  estimate  of  (3)  to  be  of  the  form 

T1  ■  “  Ki0  *  3  Kn  ’’‘i:*11!]*5  Ki4-  1  ■  1 . N  (5) 

where  K  .  is  K,.  of  the  next  integration  step.  A  reasonable  truncation 
error  eStimatewill  be  achieved  if  (5)  reflects  the  magnitude  of  the 
fourth-order  term  of  the  Taylor  series  expansions  of  (1) .  This  is  accom¬ 
plished  by  expanding  (5)  in  a  Taylor  series  and  requiring  the  first  three 
terms  to  be  tero. 

Similar  to  the  procedure  in  Appendix  A,  one  considers  the  case  of  one  de¬ 
pendent  variable;  thus  becomes 

Ki4  "  K4  "  hf(XQ  +  h>  y0  ♦  a  K0  ♦  b  K  x  ♦  c  d  Kj)  (6) 


The  error  control  is  obtained  from  a  Taylor  series  expansion  of  (6) .  For 
simplicity  of  notation,  define  operators  and  by 


k  +  <a  +  b  *  c  *  d)fn 


a 

3y 


and 


41 


h  -k  *  (4  s  ♦  '>  • c  *  -i  k3)  ~ 


=!.»,•  [b(Kj  -  hf0)  •  c(K,  -  hf0)  ♦  d(Kj  -  hfn)]  ly 
In  terms  of  operators  D ^ ,  L)., ,  and  D,  defined  in  Appendix  A, 


k4  =  h[f  ♦  d41  f 


41 


i)3 

U41 


f  *-rrf  *  ■■■) 


X  =  X, 


■  h{f  ♦  h  I),  f  ♦  L-  n/  f  ♦  h2  f  [bD.f  ♦  c  D,  f  ♦  d  D,f] 
3  432*4  y  1  2  0 

♦  ~  n,5  f  *  Vf  (bD,2  f  ♦  c  I),2  f  ♦  d  D  2  f] 

j !  4  2  vl  1  2  3 

♦  h^  f  “  [crD^f  ♦  d(sDjf  +  tl»,f)] 

• h3  c m  *  Jn3f]}x  . 


(7) 


Substituting  (A9)*  and  (7)  into  (5)  and  equating  the  first  three  terms 
of  the  Taylor  series  expansion  of  (5)  to  zero  yields  the  following  four 
algebraic  equations: 

a  +  8  +  Y  +  |5  +  e  =  0 

8m  ♦  yn  ♦  6p  +  c  -  0 

2  2  2 
8m  ♦  yn  +  6p  +  c  =  0 

yrm  ♦  <5(sm  +  tn)  ♦  e(hm  ♦  cn  +  dp)  =  0  (8) 

Since  the  four  equations  above  involve  five  new  parameters  one  may  solve 
for  them  in  terms  of  one  of  the  new  parameters.  The  general  solution  of 
(8)  in  terms  of  the  parameter  e  and  the  previously  selected  independent 
parameters  is  given  in  Table  AI  of  Appendix  A. 

To  gain  more  information  about  the  error  estimate,  it  is  required  that 
the  fourth-order  term  of  the  Taylor  series  expansion  of  (S)  match  as 
closely  as  possible  the  corresponding  fourth-order  tern  of  (AS)  .  For 
this  purpose,  fitting  coefficients  are  defined  which  express  the  difference 
between  the  coefficients  of  the  corresponding  fourth-order  terms  mentioned 
above.  By  using  (7),  (A9) ,  (A10),  (5),  and  (A5)  these  fitting  coefficients 
are  given  by 

*  1  K\  ( ♦  yv?  ♦  6p^  ♦  e)  -  1/24 

B  yrmn  +  6p  (sm  ♦  tn)  ♦  c  (bm  +  cn  ♦  dp)  -  1/8  (9) 

F^  =  1/2  [yrm"  ♦  6(sm"  ♦  tn^)  ♦  e  (hi:i“  ♦  dp^) ]  -  1/24 

F^  =  6trn  ♦  c[crn  +  d(sn  ♦  tn)  ]  -  1/24. 

The  selected  value  of  e  should  be  one  which  minimizes  these  fitting 
coefficients.  Table  AI  includes  the  expression  of  these  coefficients  in 
terms  of  the  previously  chosen  independent  parameters. 

From  the  Case  II  methods  in  Table  AI,  the  choice  of  t  =  1  yields  the 
classical  formulas  of  Range  while  t  =  1  ♦  /l/2  yields  the  Runge-kutta-Gill 
method . 

Due  to  the  simplicity  of  the  truncation  error  estimate  in  Case  II  ~'.ethods 
and  from  the  form  of  the  fitting  coefficient,  F^,  the  choice  of  Case  II 
methods  with  e  *  1  seems  to  be  advantageous. 

If  (1)  is  a  function  of  the  independent  variable  only,  the  Case  II  methods 
reduce  to  Simpson's  rule  and  the  corresponding  truncation  error  estimate 
clearly  fails. 


(Ai)  refers  to  (1)  in  Appendix  A,  i  ■  S, 6,7,8,9,10. 


2 .  Control  of  Step  Size 


To  control  the  error  in  the  numerical  computation  adequately  it  is  de¬ 
sirable  to  have  some  criterion  by  which  a  decision  may  be  made  to  alter 
the  current  integration  step  size.  For  a  given  numerical  integration 
formula  the  control  of  the  step  size  is  based  on  comparisons  of  its 
associated  truncation  error  estimates  wi:h  the  corresponding  error  limits 
for  each  variable.  Therefore,  by  keeping  the  error  limits  constant,  the 
truncation  error  will  remain  below  a  fixed  bound,  while  using  a  fixed 
percentage  of  the  current  magnitude  of  each  variable  for  the  error  limits 
keeps  the  relative  error  below  „  fixed  bound. 

In  using  the  fixed  percentage  method,  however,  oscillating  variables 
cause  the  calculations  to  bog  down  when  these  variables  cross  zero.  Thus 
a  weighted  average  of  the  previous  magnitudes  of  the  variables  could  be 
used  as  a  basis  for  the  error  limits. 


The  truncation  error  estimate,  T.,  can  be  used  to  circumvent  the  usual 

Richardson  method  of  step  size  control.  This  is  accomplished  by  assuming 

that  the  decision  as  to  whether  or  not  to  increase  the  current  step  size 

is  based  on  a  prediction  of  what  the  truncation  error  estimates  would  be 

if  a  larger  step  size  were  used.  Therefore,  let  h*  be  the  larger  step 

size  and  let  T, *  be  the  truncation  error  estimate  using  h* .  By  assuming  that 

tne  truncalivr.  error  estimate,  T.,  is  of  order  k  and  that  the  kth  term  of 

1 

the  Taylor  series  expansion  of  T.  dominates  the  expansion,  one  may  approx¬ 
imate  T\*  by  1 


T.* 

l 


h\k 


T. 


In  practice  it  is  desirable  to  overestimate  the  factor  T.*  since  this 
tends  to  reduce  time-consuming  premature  increases  in  step  size. 


For  a  fourth-order  method  with  h*  =  2h,  T.*  *  32T.,  let  U.  be  the  pre¬ 
determined  bound  for  the  truncation  erior1estinat^  T. .  If  T.  >  U.,  y. 
yi(x^  ♦!>.)  is  unacceptable .  Therefore  the  step  size  h1should  be  reduced 
and  y.  (x„  ♦h)  computed  again  using  the  smaller  value  of  h.  On  the 
other  nana,  T.  <  U.  implies  the  integration  is  acceptable.  If  in  addition, 
T.*  <  U.  ,  the' trunbat ion  error  generated  by  h*  should  be  within  an  accept¬ 
able  tolerance;  thus,  h*  is  used  as  the  integration  interval  in  the  next 
integration . 


Herein  lies  the  true  potential  of  the  above  technique  -  no  additional 
functional  evaluations  were  required  beyond  those  of  the  numerical  in¬ 
tegration  method  itself. 
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APPLICATIONS  AND  CONCLUSIONS 


The  Runge-Kutta-Gill  numerical  integration  scheme  was  incorporated  into 
the  six-degree-of-freedon  mathematical  ballistic  models  of  Walter  (4) 
and  Duncan  and  Ensey  (5).  The  coefficient  c  »  1  was  chosen  for  us*  m 
the  estimate  of  truncation  error,  several  simulated  trajectories  were 
then  computed  for  each  of  the  following  unguided  rockets:  the  Aerobee 
350  (a  high-altitude  research  rocket),  the  BlfTS  (Ballistic  Missile 
Target  System),  and  the  Athena  (a  reentry  research  rocket).  The  re¬ 
sults  of  these  simulations  were  compared  to  similar  simulations  using 
the  Runge-Kutta-Gill  integration  scheme  but  with  the  extrapolation-to- 
the-zero-grid  error  control  criterion. 

The  differences  for  each  of  the  comparable  simulations  at  burnout,  at 
peak,  and  at  impact  were  considered  negligible  (less  than  .1^).  How¬ 
ever,  there  was  one  notable  difference  -  the  mean  processing  time  for 
a  single  trajectory  simulation  was  reduced  by  46t  by  incorporating  the 
error  technique  as  developed  herein  into  the  aforementioned  ballistic 
models . 

Since  these  ballistic  models  are  to  be  used  in  real-time  support  of  un¬ 
guided  rocket  firings  where  time  is  so  essential,  there  is  a  clear  ad¬ 
vantage  in  using  the  error  check  as  developed  herein  to  replace  the 
currently  used  extrapolation-to-the-zero-grid  method 


ro 
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APPENDIX  A 


THE  GENERAL  RUNGE-KUTTA  METHOD 


Suppose  one  has  a  system  of  first-order  ordinary  differential  equations 
of  the  form 

dyt 

“dx  =  yi  =  (x,  yt,  ....  yN)  i  =  1,  ....  N  (A] 

with  initial  conditions  y.  (x^)  =  y._,  i  =  1,  N.  The  Runge-Kutta 

fourth-order  method  is  an  algorithrnaesigned  to  approximate  the  Taylor 
series  expansions  0 

(xo  * h)  *  >'io  * h  >to  *  rryio"  *  (a: 

of  (Al)  and  has  the  general  form 


whe*  e 


yi 

(X0  *  h)  = 

:  yi0 

+  aKiO 

♦ 

bK.,  ♦  c  K. _  ■< 
ll  i2 

•  dK.  _ 
i3 

(A3) 

Ki0 

= 

hf.(x0. 

yl0’ 

•  •  •,  : 

W  1  =  1* 

....  N 

Kil 

= 

hf.<*0 

+  lh. 

yio  + 

n 

K10*  •••’  yN0 

+  n  W 

(A4) 

Ki7 

= 

hfi‘x0 

+  nh. 

yio  + 

q 

K10  +  r  Kll*  ' 

....  yN0  - 

q  kno 

h  r  kni} 

Ki3 

= 

hf.(x0 

4  ph. 

yio  + 

V 

K10  +  5  Kn  + 

1  K12,  •’ 

To  reduce  somewhat  the  unwieldy  nature  of  the  notation,  one  requires  that 

3  c 

(Al)  be  a  function  of  one  dependent  variable,  and  defines  D  ■  —  +  f  — 

and  f  =  The  Taylor  series  solutions  of  Al  can  now  be  written  as 

y  3y  2  3 

y(x0  ♦  h)  =  y0  ♦  [hf  ♦  ~  Df  *  ^  (D2f  ♦  fy  Df)  (A5) 

4 

*  Tf  +  fv  °2f  +  V  Df  +  3[1f  |,fv)  +  -  r 


while  the  fourth-order  Runge-Kutta  method  tabes  on  the  following  form: 


y(x0  ♦  h)  «  y0  ♦  a  K0  ♦  bKj  +  c  K,  +  d  K, 


where 


Ko  =  hf(V  yO> 
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Kj  =  hf(x0  ♦  lh,  yQ  ♦  m  Kq) 

K2  *  hf(xQ  ♦  nh,  yQ  +  q  5C0  >  rKj) 

K3  *  hf(x0  *  Ph>  yQ  +  v  K0  +  S  K1  *  tK25 


(A7) 


D1 

°2 
D  =  p 


series 

expansion 

1  1. 

3x 

4  mf  ~ 

0  3y 

3 

n  3x 

4  (q  4  r) 

3 

P  3x 

+  (v  4  s 

"  JL 

C0  3y 

h 


(A8) 


where  fQ  =  f(xQ).  Using  A8,  the  expansions  of  K^,  and  are 

2  ,3 


K1  •  hf0 

2  3  4 

K2=h[f+hD.f^  D  3  f  ♦  Jr  D*  f  -  •  •  -]x  , 


K,  =  h  {f  +  h  D.f  ♦  ^  D  2  f  *  ~ 

•J  4m  <C  •  C  J  . 


3  V4 

°2  f  +  41  n2  f  +  ••* 


0 


♦h 


T[f  D.f  ♦  2..-  f  D.  f  4  h  D.  f  D,f  ♦  —  f  D,  f 
1  y  1  2 !  y  1  1  2y  .>!  y  1 


*  T  Ul  f  D2fy  *  5T  r  <yy  <D1  ^  *  5T  V  ^  V  !x 


K4  =  h  {f  ♦  h  D.f  ♦  h2  f  [sCDjf  ♦  yDj2  f  ♦  —  Dj3  f  +  ...)  (A9) 

2  2 

♦  t(D2  f  ♦  hr  fy  Dxf  ♦  j  D22  f  +  r  fy  Dj  2  f  4  h2  Dj  f  D2  fy  ♦  jjyD^f 

♦  jr  (^2  D32  f  *  2h3  D3  f  [s(Dj  f  +  y  Dj  2  f  4  . . .)  4  t(D2f  +  jDj2  f 


y 

4  hr  f  D.f  4  ...)]  4  h4  f  [s2(D,f)2  4  2  st  D.f  D.  f 
y  1  yyL  M  '  12 


4  t2  (D2  f)2  4  ...])  4  L-  (h3  D33  f  4  3h4  D/  f  [s  Djf  4  t  02  f  4  ...]) 


iTth4D34  f.  ...)  ♦ 
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By  requiring  that  the  Runge-Kutta  method  be  independent  of  the  choice 
of  the  function  f,  one  obtains  from  (A8) 

l=m,  n=q+r,  p=v+s+t 

i.e.  =  mD,  D2  =  nD,  *  pD,  (A10) 

Substitution  of  (A10)  into  (A9) ,  and  then  into  (A6) ,  and  equating  the 
resulting  set  of  equations  term-by-term  to  the  Taylor  series  expansion 
(AS)  yields  the  following  set  of  eight  algebraic  equations  in  ten  un¬ 
knowns  . 


a  +  b  +  c  +  d  =  l 

bm  +  cn  +  dp  -  1/2 

2  2  2 
bm  ♦cn  +  dp  =  1/3 

bm3  ♦  cn3  +  dp3  =  1/4 


crm  +  d(sm  +  tn)  =  1/6 

2  2  2 

crm  +  d(sn  +  tn  )  *  1/2 

crmn  +  dp(sm  +  tn)  =  1/8 

dtrm  =  1/24  (All) 


Table  AI  lists  the  general  solution  of  the  equations  in  (All)  in  terms  of 
the  parameters  n  and  n  and  several  special  cases  of  the  solution  which 
result  from  certain  sets  of  equations  in  (All)  being  linearly  dependent 
for  particular  values  of  m  and  n. 


» 
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TABLU  AI 

SOLUTIONS  OF  LODATIONS 


Case 

t 

ii 

III 

IV 

ifrnmrM 

m 

m 

1/2 

l 

1/2 

n 

n 

1/2 

1/2 

0 

P 

1 

1 

1 

1 

n(m-n) 

1 

1 

1 

r 

It 

F 

It 

fl-m) (m-4n2+5n-2) 

1  4. 

t 

3 

s 

2m(n-m)  C6mn-4m-4n-*-3^ 

1  -  t 

"4 

2 

(I-2m) (1-m)  (l-n) 

t 

n(n-m) (6mn-4m-4n+$3 

t 

a 

6mn-2m-2n+l 

1 

1 

l-'t 

12mn 

F 

F 

“T 

K 

2n-l 

2-t 

t-2 

2 

12m(n-m) (l-m) 

3 

6t 

T 

2m- 1 

t 

2 

t 

c 

lTn(m-n)  (l^T 

T 

I 

F 

6mn-4m-4n+3 

l 

1 

l 

d 

12(l”m)  (1-n) 

6 

3t 

6 

a 

e(2m-l)(2n-l) 

2mn 

0 

0 

-et 

B 

eC2m-l)(2n-l) 

0 

e  C  2-t J 
t 

0 

2mCn-m)  (l-m) 

u 

e(2m-l) (2n-l) 

0 

o 

Et 

Y 

2n(n-m)  (t-n) 

e(6mn-4m-4n+3) 

2c 

6 

'  '  '2fWa-n3 - 

“  £ 

t 

•  £ 

£ 

e 

E 
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